Load and visualize fits for NEON NO3 models

Author

Christa Torrens

Intro

This document uses tidybayes to extract and visualize model fit params for the NEON data/ autotrophic uptake project

First created 2/25/2025 by Christa Torrens, with ongoing amendments

Load needed libraries and model fits

Show the code
# Load libraries
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstanarm)
library(here)
library(plotly)
library(sf)

Plotting functions

This keeps all of the core plotting code in 1 place, and allows site-level customization (mostly dataframes and titles)

Show the code
##### Model fit - Npred + Nobs over time

plot_Nmodelfit_ts <- function(data, start_jday, end_jday, plot_title) {
  plot <- data %>%
    filter(model_jday >= start_jday & model_jday <= end_jday) %>%
    ggplot(aes(x = mod_hours)) +
    geom_point(aes(y = N_conc)) + 
    geom_line(aes(y = N_pred), col = 'red') +
    labs(
      x = "Time (h)",
      y = expression("N" ~ (mmol ~ m^-3))
    ) +
    ggtitle(plot_title) +
    facet_wrap(~yr_jday) +
    theme_bw()
}

plot_Nmodelfit_tsCI <- function(data, start_jday, end_jday, plot_title) {
  plot <- data %>%
    filter(model_jday >= start_jday & model_jday <= end_jday) %>%
    ggplot(aes(x = mod_hours)) +
    geom_ribbon(aes(ymin=N_pred.lower, ymax=N_pred.upper), alpha = 0.7, fill = "slategray") +
    geom_point(aes(y = N_conc)) + 
    geom_line(aes(y = N_pred), col = 'red') +
    labs(
      x = "Time (h)",
      y = expression("N" ~ (mmol ~ m^-3))
    ) +
    ggtitle(plot_title) +
    facet_wrap(~yr_jday) +
    theme_bw()
}

##### Model fit - Npred v Nobs w. 1:1 line

plot_Npred_v_Nobs <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=N_conc, y=N_pred)) +
  geom_point() + 
  labs( x=expression("measured N"~(mmol~m^-3)), # OR (mu*mol~L^-1)
        y=expression("predicted N"~(mmol~m^-3))  # ~ = a space, * = no space
      ) + 
  ggtitle(plot_title) +
  geom_abline(intercept = 0, slope = 1, color = "red") + 
  theme_bw()
}

##### Equifinality checks
# added error bars to all equifinality checks, 4/16/25

## K vs U

plot_KvU <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(y=K, x=U)) +
   geom_errorbar(aes(ymin = K.lower, ymax = K.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
   labs(y=expression("modeled K"~(d^-1)),
       x=expression("modeled U"~(mmol~m^-2~d^-1))
        ) + 
   ggtitle(plot_title) + 
  # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
   theme_bw()

}
  

## K vs N_e

plot_KvNe <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(y=K, x=N_e)) +
   geom_errorbar(aes(ymin = K.lower, ymax = K.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
   labs(y=expression("modeled K"~(d^-1)),
        x=expression("modeled N_e"~(mmol~m^-3)) 
        ) + 
   ggtitle(plot_title) + 
   # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
   theme_bw()
}


## U vs N_e
# Includes CI bars because this relationship tends to look a bit linear - checking whether that holds even w ci values in place

plot_UvNe <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(y=U, x=N_e)) +
   geom_errorbar(aes(ymin = U.lower, ymax = U.upper), 
                width = 0.2, color='grey70') + # CI bars
   geom_point() + 
   labs(y=expression("modeled U"~(mmol~m^-2~d^-1)),
        x=expression("modeled N_e"~(mmol~m^-3)) 
        ) + 
   ggtitle(plot_title) + 
   theme_bw()
}


##### K over time

plot_K_ts <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=d, y=K)) +
   geom_errorbar(aes(ymin = K.lower, ymax = K.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
  # geom_point(y = sumlight.real, color = 'gold') +
   labs( x= "Index", 
         y=expression("modeled K"~(d^-1)) # ~ = a space, * = no space
      ) + 
   ggtitle(plot_title) +
   theme_bw()
}

  
##### U over time

plot_U_ts <- function(data, plot_title) {
  ggplot(data=data, aes(x=d, y=U)) +
   geom_errorbar(aes(ymin = U.lower, ymax = U.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
  # geom_point(y = sumlight.real_wlou, color = 'gold') +
  # ADD IN HIGH AND LOW CIs
   labs( x= "Index", 
        y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
      ) + 
   ggtitle(plot_title) +
   theme_bw()
}


##### U vs. real sumlight

plot_UvLight <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=sumlightReal, y=U)) +
  geom_point() + 
  labs( x=expression("measured daily light"~(w~m^-2)), 
        y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
      ) + 
  scale_x_log10() + scale_y_log10() + 
  ggtitle(plot_title) +
  theme_bw()
}


##### U + daily (or hourly?) sumlight over time



##### Autotrophic U vs total U

plot_autoUperc <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=d, y=U_auto_perc)) +
  geom_point() + 
  xlab("Index") + ylab("Percent autotrophic uptake") + 
  ggtitle(plot_title) + 
  theme_bw()
}

Load model fits and data by site (stored as RDS files), then extract parameters

BIGC

Outputs and visualizations for Upper Big Creek, Fresno, CA

Load model fits and data
Show the code
########## Load model fit ______________________________________________________

bigc.fit <- readRDS(here("N_uptake_NEON/data/model_fits/bigc.fit.rds"))

# if working from model run
# bigc.fit <- fit.bigc


########## Load model data _____________________________________________________

bigc.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/bigc.data.rds"))


########## Load full daily and hourly datasets _________________________________
# these include datatime info, and the model data does not)

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/bigc_clean.csv")
# bigc.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="US/Pacific"), 
#          model_datetime = local_datetime - hours(4))  


# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/bigc_hourly_clean.csv")
bigc.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Pacific"), 
         model_datetime = local_datetime - hours(4))  
Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a BIGC tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd   2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.05    0.00 0.00   0.05  0.05  0.05  0.05  0.06  2811 1.00
sigma_U  0.16    0.00 0.04   0.09  0.13  0.16  0.18  0.23   129 1.02
b0      -9.15    0.03 0.74 -10.69 -9.62 -9.14 -8.67 -7.77   667 1.01
b1       0.82    0.00 0.08   0.66  0.77  0.82  0.88  1.00   683 1.01

Samples were drawn using NUTS(diag_e) at Tue Apr  1 09:45:41 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
Create visualizations
Show the code
########## N Model fit to data __________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_bigc <- plot_Nmodelfit_ts(data=N_output_bigc.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Upper Big Creek, CA")

# quartz(width=6,height=6)
N_Npred_ts_bigc

Show the code
### w. credible interval ribbon
N_Npred_tsCI_bigc <- plot_Nmodelfit_tsCI(data=N_output_bigc.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Upper Big Creek, CA")

# quartz(width=6,height=6)
N_Npred_tsCI_bigc

Show the code
######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_bigc <- plot_Npred_v_Nobs(data=N_output_bigc.df, plot_title="Predicted N vs. measured N, Upper Big Creek, CA")

# quartz(width=6, height=6)
Npred_v_N_bigc

Show the code
########## Equifinality checks _________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_bigc <- plot_KvU(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_bigc

Show the code
####### Check equifinality!! K vs N_e  #######  

KvNe_bigc <- plot_KvNe(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_bigc

Show the code
####### Check equifinality!! U vs N_e  #######  

UvNe_bigc <- plot_UvNe(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_bigc

Show the code
########## Other parameter visualizations _______________________________________

#######  K over time  #######  

K_v_time_bigc <- plot_K_ts(data=daily_output_bigc.df, plot_title = "modeled K over time, Upper Big Creek, CA")

# quartz(width=6.5, height=6)
K_v_time_bigc

Show the code
##### U over time  #######  

U_time_bigc <- plot_U_ts(data=daily_output_bigc.df, plot_title="Diel nitrate uptake (modeled), Upper Big Creek, CA")

# quartz(width=6, height=6)
U_time_bigc

Show the code
###### U vs sumlight  #######  

U_vs_light_bigc <- plot_UvLight(data=daily_output_bigc.df, plot_title = "Scatterplot of NO3 uptake and daily light: Upper Big Creek, CA")

# quartz(width=6.5, height=6)
U_vs_light_bigc

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

Show the code
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_bigc.df <- daily_output_bigc.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_bigc.df$U_auto_perc) #11.7% on average

Uauto_perc_bigc <- plot_autoUperc(data=daily_output_bigc.df, plot_title = "Daily  autotrophic uptake (as % of total uptake), Upper Big Creek, CA")

# quartz(width=7.5, height=6)
Uauto_perc_bigc 

Show the code
#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_bigc)
Save updated model output datasets
Show the code
# path = here("N_uptake_NEON/data/model_output/bigc_output_daily")
# write_csv(daily_output_bigc.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/bigc_output_hourly")
# write_csv(N_output_bigc.df, file=path.h)

CARI

Outputs and visualizations for Caribou Creek, Fairbanks North Star, AK

Load model fits and data
Show the code
########## Load model fit ______________________________________________________

cari.fit <- readRDS(here("N_uptake_NEON/data/model_fits/cari.fit.rds"))

# if working from model run
# cari.fit <- fit.cari


########## Load model data _____________________________________________________

cari.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cari.data.rds"))


########## Load full hourly datasets __________________________________
# these include datatime info, and the model data does not)

# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/cari_hourly_clean.csv")
cari.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Alaska"), 
         model_datetime = local_datetime - hours(2))# 2-hr offset because it starts getting light at 02:45 in midsummer...
Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a CARI tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=3000; warmup=1500; thin=1; 
post-warmup draws per chain=1500, total post-warmup draws=6000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.10       0 0.00  0.10  0.10  0.10  0.11  0.11 10850    1
sigma_U  0.28       0 0.02  0.23  0.26  0.28  0.29  0.32  1579    1
b0      -2.26       0 0.29 -2.85 -2.45 -2.25 -2.06 -1.69  3517    1
b1       0.19       0 0.03  0.13  0.17  0.19  0.21  0.25  3516    1

Samples were drawn using NUTS(diag_e) at Tue Apr 15 15:00:35 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
Create visualizations - Caribou Creek, AK
Show the code
########## N Model fit to data _________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_cari <- plot_Nmodelfit_ts(data=N_output_cari.df, start_jday = 145, end_jday = 155, plot_title = "NEON N data + model predictions: Caribou Creek, AK")

# quartz(width=6,height=6)
N_Npred_ts_cari

Show the code
### w. credible interval ribbon
N_Npred_tsCI_cari <- plot_Nmodelfit_tsCI(data=N_output_cari.df, start_jday = 145, end_jday = 155, plot_title = "NEON N data + model predictions: Caribou Creek, AK")


# quartz(width=6,height=6)
N_Npred_tsCI_cari

Show the code
######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_cari <- plot_Npred_v_Nobs(data=N_output_cari.df, plot_title="Predicted N vs. measured N, Caribou Creek, AK")

# quartz(width=6, height=6)
Npred_v_N_cari

Show the code
########## Equifinality checks ________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_cari <- plot_KvU(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_cari

Show the code
####### Check equifinality!! K vs N_e  #######  

KvNe_cari <- plot_KvNe(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_cari

Show the code
####### Check equifinality!! U vs N_e  #######  

UvNe_cari <- plot_UvNe(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_cari

Show the code
########## Other parameter visualizations ______________________________________

#######  K over time  #######  

K_v_time_cari <- plot_K_ts(data=daily_output_cari.df, plot_title = "modeled K over time, Caribou Creek, AK")

# quartz(width=6.5, height=6)
K_v_time_cari

Show the code
##### U over time  #######  

U_time_cari <- plot_U_ts(data=daily_output_cari.df, plot_title="Diel nitrate uptake (modeled), Caribou Creek, AK")

# quartz(width=6, height=6)
U_time_cari

Show the code
###### U vs sumlight  #######  

U_vs_light_cari <- plot_UvLight(data=daily_output_cari.df, plot_title = "Scatterplot of NO3 uptake and daily light: Caribou Creek, AK")

# quartz(width=6.5, height=6)
U_vs_light_cari

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

Show the code
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_cari.df <- daily_output_cari.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_cari.df$U_auto_perc) #4.9% on average

Uauto_perc_cari <- plot_autoUperc(data=daily_output_cari.df, plot_title = "Daily  autotrophic uptake (as % of total uptake), Caribou Creek, AK")

quartz(width=7.5, height=6)
Uauto_perc_cari 

#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_cari)
Save updated model output datasets
Show the code
# path = here("N_uptake_NEON/data/model_output/cari_output_daily")
# write_csv(daily_output_cari.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/cari_output_hourly")
# write_csv(N_output_cari.df, file=path.h)

CUPE

Outputs and visualizations for Rio Cupeyes, San German Municipio, PR

Load model fits and data
Show the code
########## Load model fit ______________________________________________________

cupe.fit <- readRDS(here("N_uptake_NEON/data/model_fits/cupe.fit.rds"))

# if working from model run
# cupe.fit <- fit.cupe


########## Load model data _____________________________________________________

cupe.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cupe.data.rds"))

########## Load full daily and hourly datasets __________________________________
# these include datatime info, and the model data does not)

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/cupe_clean.csv")
# cupe.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="America/Puerto_Rico"), 
#          model_datetime = local_datetime - hours(4))  


# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/cupe_hourly_clean.csv")
cupe.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="America/Puerto_Rico"), 
         model_datetime = local_datetime - hours(4))  
Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a CUPE tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.15    0.00 0.00  0.14  0.15  0.15  0.15  0.15  5813    1
sigma_U  0.24    0.00 0.01  0.22  0.24  0.24  0.25  0.26  2620    1
b0      -5.25    0.01 0.44 -6.08 -5.56 -5.25 -4.96 -4.37  2832    1
b1       0.56    0.00 0.05  0.46  0.53  0.56  0.60  0.66  2841    1

Samples were drawn using NUTS(diag_e) at Tue Apr  1 15:52:19 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
Create visualizations - Rio Cupeyes, San German Municipio, PR
Show the code
########## N Model fit to data __________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_cupe <- plot_Nmodelfit_ts(data=N_output_cupe.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Rio Cupeyes, PR")

# quartz(width=6,height=6)
N_Npred_ts_cupe

Show the code
### w. credible interval ribbon
N_Npred_tsCI_cupe <- plot_Nmodelfit_tsCI(data=N_output_cupe.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Rio Cupeyes, PR")

# quartz(width=6,height=6)
N_Npred_tsCI_cupe

Show the code
######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_cupe <- plot_Npred_v_Nobs(data=N_output_cupe.df, plot_title="Predicted N vs. measured N, Rio Cupeyes, PR")

# quartz(width=6, height=6)
Npred_v_N_cupe

Show the code
########## Equifinality checks _________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_cupe <- plot_KvU(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_cupe

Show the code
####### Check equifinality!! K vs N_e  #######  

KvNe_cupe <- plot_KvNe(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_cupe

Show the code
####### Check equifinality!! U vs N_e  #######  

UvNe_cupe <- plot_UvNe(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_cupe

Show the code
########## Other parameter visualizations ______________________________________

#######  K over time  #######  

K_v_time_cupe <- plot_K_ts(data=daily_output_cupe.df, plot_title = "modeled K over time, Rio Cupeyes, PR")

# quartz(width=6.5, height=6)
K_v_time_cupe

Show the code
##### U over time  #######  

U_time_cupe <- plot_U_ts(data=daily_output_cupe.df, plot_title="Diel nitrate uptake (modeled), Rio Cupeyes, PR")

# quartz(width=6, height=6)
U_time_cupe

Show the code
###### U vs sumlight  #######  

U_vs_light_cupe <- plot_UvLight(data=daily_output_cupe.df, plot_title = "Scatterplot of NO3 uptake and daily light: Rio Cupeyes, PR")

# quartz(width=6.5, height=6)
U_vs_light_cupe

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

Show the code
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K * z
# TOTAL UPTAKE = U + (Equilibrium nitrate * K * z)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_cupe.df <- daily_output_cupe.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_cupe.df$U_auto_perc) #4.7% on average

Uauto_perc_cupe <- plot_autoUperc(data=daily_output_cupe.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Rio Cupeyes, PR")

# quartz(width=7.5, height=6)
Uauto_perc_cupe 

Show the code
#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_cupe)
Save updated model output datasets
Show the code
# path = here("N_uptake_NEON/data/model_output/cupe_output_daily")
# write_csv(daily_output_cupe.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/cupe_output_hourly")
# write_csv(N_output_cupe.df, file=path.h)

PRIN

Outputs and visualizations for Pringle Creek, Wise County, TX

Load model fits and data
Show the code
########## Load model fit ______________________________________________________

prin.fit <- readRDS(here("N_uptake_NEON/data/model_fits/prin.fit.rds"))

# if working from model run
# prin.fit <- fit.prin


########## Load model data _____________________________________________________

prin.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/prin.data.rds"))

########## Load full daily and hourly datasets __________________________________
# these include datatime info, and the model data does not)

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/prin_clean.csv")
# prin.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="US/Central"), 
#          model_datetime = local_datetime - hours(4))  


# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/prin_hourly_clean.csv")
prin.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Central"), 
         model_datetime = local_datetime - hours(4))  
Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a PRIN tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.31    0.00 0.00  0.30  0.30  0.31  0.31  0.32  6493    1
sigma_U  0.36    0.00 0.04  0.28  0.33  0.35  0.38  0.44  1104    1
b0      -6.48    0.03 1.26 -9.11 -7.28 -6.48 -5.65 -4.03  1608    1
b1       0.70    0.00 0.15  0.42  0.61  0.70  0.80  1.00  1621    1

Samples were drawn using NUTS(diag_e) at Wed Apr 16 00:03:30 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
Create visualizations - Pringle Creek, Wise County, TX
Show the code
########## N Model fit to data __________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_prin <- plot_Nmodelfit_ts(data=N_output_prin.df, start_jday = 90, end_jday = 120, plot_title = "NEON N data + model predictions: Pringle Creek, TX")

# quartz(width=6,height=6)
N_Npred_ts_prin

Show the code
### w. credible interval ribbon
N_Npred_tsCI_prin <- plot_Nmodelfit_tsCI(data=N_output_prin.df, start_jday = 90, end_jday = 120, plot_title = "NEON N data + model predictions: Pringle Creek, TX")

# quartz(width=6,height=6)
N_Npred_tsCI_prin

Show the code
######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_prin <- plot_Npred_v_Nobs(data=N_output_prin.df, plot_title="Predicted N vs. measured N, Pringle Creek, TX")

# quartz(width=6, height=6)
Npred_v_N_prin

Show the code
########## Equifinality checks _________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_prin <- plot_KvU(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_prin

Show the code
####### Check equifinality!! K vs N_e  #######  

KvNe_prin <- plot_KvNe(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_prin

Show the code
####### Check equifinality!! U vs N_e  #######  

UvNe_prin <- plot_UvNe(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_prin

Show the code
########## Other parameter visualizations ______________________________________

#######  K over time  #######  

K_v_time_prin <- plot_K_ts(data=daily_output_prin.df, plot_title = "modeled K over time, Pringle Creek, TX")

# quartz(width=6.5, height=6)
K_v_time_prin

Show the code
##### U over time  #######  

U_time_prin <- plot_U_ts(data=daily_output_prin.df, plot_title="Diel nitrate uptake (modeled), Pringle Creek, TX")

# quartz(width=6, height=6)
U_time_prin

Show the code
###### U vs sumlight  #######  

U_vs_light_prin <- plot_UvLight(data=daily_output_prin.df, plot_title = "Scatterplot of NO3 uptake and daily light: Pringle Creek, TX")

# quartz(width=6.5, height=6)
U_vs_light_prin

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

Show the code
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K * z
# TOTAL UPTAKE = U + (Equilibrium nitrate * K * z)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_prin.df <- daily_output_prin.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_prin.df$U_auto_perc) #12.74% on average

Uauto_perc_prin <- plot_autoUperc(data=daily_output_prin.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Pringle Creek, TX")

# quartz(width=7.5, height=6)
Uauto_perc_prin 

Show the code
#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_prin)
Save updated model output datasets
Show the code
# path = here("N_uptake_NEON/data/model_output/prin_output_daily")
# write_csv(daily_output_prin.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/prin_output_hourly")
# write_csv(N_output_prin.df, file=path.h)

WLOU

Ouptuts and visualizations for West St Louis Creek, Grand, CO

Load model fits and data
Show the code
########## Load model fit ______________________________________________________

wlou.fit <- readRDS(here("N_uptake_NEON/data/model_fits/wlou.fit.rds"))

# if working from model run
# wlou.fit <- fit.wlou


########## Load model data _____________________________________________________

# data from model

wlou.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/wlou.data.rds"))

###### Load dataset (includes datatime info, which model data does not) _______

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/wlou_clean.csv")
# wlou.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="US/Mountain"), 
#          model_datetime = local_datetime - hours(4))  
# 

# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/wlou_hourly_clean.csv")
wlou.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Mountain"), 
         model_datetime = local_datetime - hours(4))  
Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a WLOU tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.05    0.00 0.00  0.05  0.05  0.05  0.05  0.05  5201    1
sigma_U  0.41    0.00 0.03  0.35  0.39  0.41  0.43  0.47  2253    1
b0      -2.64    0.01 0.57 -3.74 -3.03 -2.64 -2.25 -1.52  3282    1
b1       0.11    0.00 0.07 -0.02  0.07  0.11  0.16  0.24  3283    1

Samples were drawn using NUTS(diag_e) at Mon Mar 10 15:19:42 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
Create visualizations - West St. Louis Creek, CO
Show the code
########## N Model fit to data ________________________________________________

######  N model fit over time  #######

# N_and_Npred_wlou <- N_output_wlou.df %>%
#   #group_by(year(model_datetime_wlou)) %>%
#   filter(model_jday >= 35 & model_jday <= 60) %>%  # to see the little boxes better...
#   ggplot(aes(x=mod_hours)) +
#   geom_point(aes(y=N_conc)) + 
#   geom_line(aes(y=N_conc_pred), col='red')+
#   labs(
#     x="Time (h)", y=expression("N"~(mmol~m^-3))
#   ) +
#   ggtitle("NEON: West St. Louis Creek N data + model predictions") +
#   facet_wrap(~yr_jday) +
#   theme_bw()

N_Npred_ts_wlou <- plot_Nmodelfit_ts(data=N_output_wlou.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: West St. Louis Creek, CO ")

#quartz(width=6,height=6)
N_Npred_ts_wlou

Show the code
### w. credible interval ribbon
N_Npred_tsCI_wlou <- plot_Nmodelfit_tsCI(data=N_output_wlou.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: West St. Louis Creek, CO")

# quartz(width=6,height=6)
N_Npred_tsCI_wlou

Show the code
######  N-pred vs N  #######

# Npred_v_N_wlou <- ggplot(data = N_output_wlou.df, aes(x=N_conc, y=N_conc_pred)) +
#   geom_point() + 
#   labs( x=expression("measured N"~(mu*mol~L^-1)), # fcn changed both to mmol~m^-3
#         y=expression("modeled N"~(mu*mol~L^-1))  # ~ = a space, * = no space
#       ) + 
#   ggtitle("Measured N vs predicted N, West St. Louis Creek, CO") +
#   geom_abline(intercept = 0, slope = 1, color = "red") + 
#   theme_bw()

Npred_v_N_wlou <- plot_Npred_v_Nobs(data=N_output_wlou.df, plot_title="Predicted N vs. measured N, West St. Louis Creek, CO")

#quartz(width=6, height=6)
Npred_v_N_wlou

Show the code
########## Equifinality checks ________________________________________________

####### Check equifinality!! K vs U  #######  

# KvU_wlou <- ggplot(data=daily_output_wlou.df, aes(y=K_mod_avg_wlou, x=U_mod_avg_wlou)) +
#   geom_point() + 
#   labs(y=expression("modeled K"~(d^-1)),
#        x=expression("modeled U"~(mmol~m^-2~d^-1))
#         ) + 
#   ggtitle("Equifinality check, West St. Louis Creek, CO: modeled K vs modeled U") + 
#   # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
#   theme_bw()


KvU_wlou <- plot_KvU(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled K vs modeled U")

  
#quartz(width=7.5, height=6)
KvU_wlou

Show the code
####### Check equifinality!! K vs N_e  #######  

# KvNe_wlou <- ggplot(data=daily_output_wlou.df, aes(y=K, x=N_e)) +
#    geom_point() + 
#    labs(y=expression("modeled K"~(d^-1)),
#         x=expression("modeled N_e"~(mmol~m^-3)) 
#         ) + 
#    ggtitle("Equifinality check, West St. Louis Creek, CO: modeled K vs modeled equilibrium N") + 
#    # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
#   theme_bw()


KvNe_wlou <- plot_KvNe(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled K vs modeled equilibrium N")

#quartz(width=7.5, height=6)
KvNe_wlou

Show the code
####### Check equifinality!! U vs N_e  #######  

# UvNe_wlou <- ggplot(data=daily_output_wlou.df, aes(y=U, x=N_e)) +
#    geom_point() + 
#    labs(y=expression("modeled U"~(mmol~m^-2~d^-1)),
#         x=expression("modeled N_e"~(mmol~m^-3)) 
#         ) + 
#    ggtitle("Equifinality check, West St. Louis Creek, CO: modeled U vs modeled equilibrium N") + 
#   theme_bw()

UvNe_wlou <- plot_UvNe(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled U vs modeled equilibrium N")
 
#quartz(width=7.5, height=6)
UvNe_wlou

Show the code
########## Other parameter visualizations ______________________________________

######## K over time  #######  

# K_v_time_wlou <- ggplot(data = daily_output_wlou.df, aes(x=d, y=K)) +
#   geom_point() + 
#   # geom_point(y = sumlight.real, color = 'gold') +
#   # ADD IN HIGH AND LOW CIs
#   # xlab("Julian day") + ylab("modeled K (day -1)") + # daily change in N concentration
#   labs( x= "Index", 
#         y=expression("modeled K"~(d^-1)) # ~ = a space, * = no space
#       ) + 
#   #ylim(0,20) +
#   ggtitle("modeled K over time, West St. Louis Creek, CO") +
#   theme_bw()


K_v_time_wlou <- plot_K_ts(data=daily_output_wlou.df, plot_title = "modeled K over time, West St. Louis Creek, CO")

#quartz(width=6.5, height=6)
K_v_time_wlou

Show the code
##### to clip: no function yet
# K_time_clip_wlou <- ggplot(data = daily_output_wlou.df, aes(x=mod_day_wlou, y=K_mod_avg_wlou)) +
#   geom_point() + 
#   # geom_point(y = sumlight.real, color = 'gold') +
#   # ADD IN HIGH AND LOW CIs
#   xlab("Julian day") + ylab("modeled K (day -1)") + # ??? UNITS?
#   ylim(0,10) +
#   ggtitle("modeled K over time, West St. Louis Creek, CO - ylim = 0,10") +
#   theme_bw()
# 
# quartz()
# K_time_clip_wlou


######## U over time  #######  

# U_time_wlou <- ggplot(data = daily_output_wlou.df, aes(x=d, y=U)) +
#   geom_point() + 
#   # geom_point(y = sumlight.real_wlou, color = 'gold') +
#   # ADD IN HIGH AND LOW CIs
#   #xlab("Julian day") + ylab("modeled U (mmol/m2/day)") + 
#   labs( x= "Index", 
#         y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
#       ) + 
#   #ylim(0,1) +
#   #ggtitle("modeled U over time, Big Creek 2019 pooled model w real light") +
#   ggtitle("Diel nitrate uptake (modeled), West St. Louis Creek, CO") +
#   theme_bw()

U_time_wlou <- plot_U_ts(data=daily_output_wlou.df, plot_title="Diel nitrate uptake (modeled), West St. Louis Creek, CO")

#quartz(width=6, height=6)
U_time_wlou

Show the code
# U_time_clip_wlou <- ggplot(data = daily_output_wlou.df, aes(x=mod_day_wlou, y=U_mod_avg_wlou)) +
#    geom_point() +
#    # geom_point(y = sumlight.real, color = 'gold') +
#    # ADD IN HIGH AND LOW CIs
#    xlab("Julian day") + ylab("modeled U (mmol/m2/day)") +
#    ylim(0,1) +
#    ggtitle("modeled U over time, West St. Louis Creek, CO - ylim = 0-1") +
#    theme_bw()
# 
#  quartz()
#  U_time_clip_wlou



######## U vs sumlight  #######  

# U_vs_light_wlou <- ggplot(data = daily_output_wlou.df, aes(x=sumlightReal, y=U)) +
#   geom_point() + 
#   #xlab("true light (satellite)") + ylab("modeled U (mmol/m2/day)") +
#   # xlab("light (satellite)") + ylab("modeled U (mmol/m2/day)") +
#   labs( x= "light (satellite)", 
#         y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
#       ) + 
#   scale_x_log10() + scale_y_log10() + 
#   ggtitle("Scatterplot of NO3 uptake and daily light: West St. Louis Creek, CO") +
#   #ylim = c(-0.2, 1) +
#   theme_bw()

U_vs_light_wlou <- plot_UvLight(data=daily_output_wlou.df, plot_title = "Scatterplot of NO3 uptake and daily light: West St. Louis Creek, CO")

# quartz(width=6.5, height=6)
U_vs_light_wlou

Compare U_total and U_auto
Show the code
# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_wlou.df <- daily_output_wlou.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)


U_auto_avg <- mean(daily_output_wlou.df$U_auto_perc) #10.3% on average

# Uauto_plot_wlou <- ggplot(data=daily_output_wlou.df, aes(x=d, y=U_auto_perc)) +
#   geom_point() + 
#   xlab("Index") + ylab("Percent autotrophic uptake") + 
#   ggtitle("Percent of autotrophic uptake each day, West St. Louis Creek, CO") + 
#   theme_bw()


Uauto_perc_wlou <- plot_autoUperc(data=daily_output_wlou.df, plot_title = "Daily  autotrophic uptake (as % of total uptake), West St. Louis Creek, CO")

# quartz(width=7.5, height=6)
Uauto_perc_wlou 

Show the code
#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_wlou)
Save updated model output datasets
Show the code
# path = here("N_uptake_NEON/data/model_output/wlou_output_daily")
# write_csv(daily_output_wlou.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/wlou_output_hourly")
# write_csv(N_output_wlou.df, file=path.h)

Assemble and save mean, median, 95% QI tibble for all sites

Combine the site tibbles into one large tibble, then save it

Show the code
########## Reload all datasets as needed  _______________________________________


########## Bind then save daily parameters (K, U, N_e)  _________________________

# daily_summary_all.df <- bind_rows()
# path <- here()
# saveRDS(daily_summary_all.df, path)

########## Bind then save hourly parameters (N_pred)  ___________________________

# hourly_summary_all.df <- bind_rows()
# path2 <- here()
# saveRDS(hourly_summary_all.df, path2)

Create whole-dataset visualizations

Show the code
########## Violin plots _______________________________________________

##### K


##### U


##### N_e